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High-Order Space-Time Methods for Conservation Laws 


H.T. Huynh 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract. Current high-order methods such as discontinuous Galerkin and/or flux reconstruction can 
provide effective discretization for the spatial derivatives. Together with a time discretization, such 
methods result in either too small a time step size in the case of an explicit scheme or a very large system 
in the case of an implicit one. To tackle these problems, two new high-order space-time schemes for 
conservation laws are introduced: the first is explicit and the second, implicit. The explicit method here, 
also called the moment scheme, achieves a CFL (Courant-Friedrichs-Lewy) condition of 1 for the case of 
one-spatial dimension regardless of the degree of the polynomial approximation. (For standard explicit 
methods, if the spatial approximation is of degree p, then the time step sizes are typically proportional to 
1/p 2 .) Fourier analyses for the one and two-dimensional cases are carried out. The property of super 
accuracy (or super convergence) is discussed. The implicit method is a simplified but optimal version of 
the discontinuous Galerkin scheme applied to time. It reduces to a collocation implicit Runge-Kutta (RK) 
method for ordinary differential equations (ODE) called Radau IIA. The explicit and implicit schemes are 
closely related since they employ the same intermediate time levels, and the former can serve as a key 
building block in an iterative procedure for the latter. A limiting technique for the piecewise linear 
scheme is also discussed. The technique can suppress oscillations near a discontinuity while preserving 
accuracy near extrema. Preliminary numerical results are shown. 

1. Introduction 

The discontinuous Galerkin (DG) method is currently among the most widely used high-order 
numerical methods for solving the compressible Navier-Stokes equations on unstructured meshes. It was 
introduced for the neutron transport equation by Reed and Hill (1973), analyzed by LaSaint and Raviart 
(1974) and developed and made popular for fluid dynamics equations by Cockbum, Shu, Bassi, Rebay, 
and others (see e.g., Bassi and Rebay 1997a,b, Cockbum, Kamiadakis, and Shu 2000, Cockburn and Shu 
2005, 2009, Shu 2012, and the references therein). Efficient DG schemes using the Lagrange polynomials 
for nodal points as basis functions known as nodal DG methods can be found in (Hesthaven and 
Warburton 2008). Alternative methods employing the differential form (as opposed to DG, which uses the 
integral form) have been proposed. Kopriva and Kolias (1996) pioneered this approach with the 
staggered-grid method on quadrilateral meshes. The method was extended to triangular meshes by Liu, 
Vinokur, and Wang (2006) and named spectral difference (SD). Another class of schemes called spectral 
volume (SV) was presented by Wang, Zhang, and Liu (2004). 

Recently, an approach to high-order accuracy with the advantage of simplicity and economy called 
flux reconstruction (FR) was introduced in (Huynh 2007, 2009a). The approach employs the differential 
form and results in several new schemes with favorable properties. In addition, by recovering DG, SD, 
and SV, it unifies existing methods. Extensions of the FR method to 2D unstructured meshes were 
presented in (Wang and Gao 2009) and (Gao and Wang 2009) and was named lifting collocation penalty 
or LCP schemes. The involved authors later combined FR and LCP and called them correction procedure 
via reconstruction or CPR methods. Extension to 3D was presented in (Haga, Gao, and Wang 2010, 201 1) 
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and (Wang, Gao, and Haga 2011). A mathematical foundation for FR was provided via energy 
stability proofs: for an SD scheme identified via FR by Jameson (2010), for a one -parameter class of FR 
schemes in ID by Vincent, Castonguay, and Jameson (2011), for 2D triangular grids by Castonguay, 
Vincent and Jameson (2012), for linear-advection-diffusion problem in ID by Castonguay et al. (2013) 
and for advection-diffusion in 2D, by Williams et al. (2011). Extension for a particularly simple FR 
scheme called g 2 to the Navier-Stokes equations was presented in (Liang, Miyaji, and Zhang 2013). 
Applications of FR/CPR methods to turbulent internal flows for turbomachinaries were elaborated in (Lu, 
Yuan, and Dawes 2012 and Lu, Liu, and Dawes 2013). Interface elements dealing with non-conforming 
polynomials together with p-adaptation for viscous flow simulations were presented in (Cagnone and 
Nadarajah 2012, Cagnone, Vermeire, and Nadarajah 2013), and an implicit large eddy simulation (ILES) 
solver was developed for FR/CPR schemes by Vermeire, Cagnone, and Nadarajah (2013). 

The high-order methods discussed can provide effective discretization for the spatial derivatives. 
Together with a time discretization, however, such methods result in either too small a time step size in 
the case of an explicit scheme or a very large system in the case of an implicit one. Efforts to improve 
time stepping have been proposed by several authors. 

Concerning explicit schemes, a space-time DG method with a favorable CFL (Courant-Friedrichs- 
Lewy) condition employing a staggered mesh was introduced in Lowrie et. al (1995). For the case of a 
semi-discrete formulation, Warburton and Hagstrom (2008) presented a staggered-mesh scheme where 
the solution is projected onto a mesh staggered from the original one and then projected back onto the 
original mesh resulting in a solution ‘smoother’ than the original. They obtain a CEL condition 
proportional to 1/p as opposed to 1/p 2 where p is the degree of the approximating polynomials. 
However, multi-dimensional extensions of these staggered-mesh approaches are both complex and costly. 
Additional efforts can be found under the heading “Taming the CEL number” of (Hesthaven and 
Warburton 2008) and the references therein. In addition to these efforts, Dumbser and Munz (2005) 
developed ADER-DG using the ADER finite-volume scheme (Arbitrary order using derivatives). Space- 
time expansion DG (or STE-DG) as a modification of ADER-DG was introduced in (Lorcher, Gassner, 
and Munz 2007 and Gassner, Lorcher, and Munz 2008), but these methods have highly restricted CEL 
limits. 

Researchers in the finite-volume community have devised several explicit methods with a CFL 
condition of 1 for the case of one spatial dimension. In fact, in 1977, Van Leer presented a series of five 
schemes for convection with this property. Among the five, scheme I is the least accurate but becomes the 
most popular and, in the literature, is generally the scheme implied when the MUSCL approach is 
referred. Concerning the more accurate schemes such as 111 and V, the problem is, as stated in (Van Leer 
and Nomura 2005), “When trying to extend these schemes beyond advection, viz., to a nonlinear 
hyperbolic system like the Euler equations, the first author ran into insuperable difficulties because the 
exact shift operator no longer applies, and he abandoned the idea”. 

For scheme III, which can be considered as a piecewise-linear space-time DG method for convection, 
this difficulty was overcome in (Huynh 2006) where the exact shift operator was replaced by (a) an 
integration by parts (Gauss theorem), (b) the observation that the space-time fluxes contain the 
information provided by the shift operator, and (c) a successive procedure of updating the moments. The 
resulting method is called the upwind moment scheme. The approach was further analyzed and applied to 
nonlinear hyperbolic equations by Suzuki and Van Leer (2007). As briefly discussed in (Huynh 2007), the 
moment scheme can be extended to arbitrary order. Such an extension was independently obtained by Lo 
and was studied in combination with Van Leer’s recovery scheme for diffusion (2011). The description of 
the extension here and that of Lo thus have considerable similarities (the differences will also be 
discussed). Extensions for scheme V can be found in (Eymann and Roe 2013). 

Concerning implicit schemes, it was shown by LaSaint and Raviart (1974) via a quadrature of left 
Radau type (Hildebrand 1987) that discretizing an ODE via discontinuous Galerkin results in an implicit 
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Runge-Kutta (RK) method. In addition, the scheme is accurate to order 2p + 1 if the solution is 
approximated by a polynomial of degree p in time. Implicit space-time DG schemes for conservation 
laws have been studied by several authors (e.g., Bar-Yoseph and Elata 1990, Van der Vegt and Van der 
Ven 2002, Pesch and Van der Vegt 2008). The common difficulty is that the implicit system is large. 
Recent effort includes the large eddy simulation solver by Vermeire, Cagnone, and Nadarajah (2013). 

Another aspect of time stepping methods concerns their performance for large scale calculations using 
massively parallel supercomputers or graphics processing units (GPU). Here, the DG or FR/CPR spatial 
discretization is highly parallelizable. It is well known that this ease of parallelization is retained with an 
explicit time stepping, e.g., RK, but not an implicit one. 

In this paper, to tackle the time stepping problem, two new high-order space-time schemes for 
conservation laws are introduced: the first is explicit and the second, implicit. The explicit method, called 
the moment scheme, achieves a CFL condition of 1 for the case of one-spatial dimension regardless of the 
degree of approximation. Compared with standard approaches, the key difference is the built-in physics of 
advection. Such a CFL condition is a significant gain compared with a typical time step size proportional 
to 1/p 2 . For the case of multiple-dimensions, however, there is a reduction in the CFL limits. The current 
extension differs from that of Lo in (a) the case of two-spatial dimensions and its Fourier analysis are 
presented and (b) additional details, observations, and the relation with the implicit time stepping are 
discussed. It will be shown via Fourier analysis that for both one and two-dimensional cases, the moment 
schemes are accurate to order 2p + 1 if the solution is approximated by polynomials of degree p. 

The implicit method here is a simplified but optimal DG scheme applied to time. It employs the FR 
approach and simplification is made possible by the fact that time always moves forward (as opposed to 
space where waves can travel in all directions). The method reduces to an implicit RK scheme using the 
right Radau points as collocation points called the Radau IIA (Lambert 1991, Hairer, Norsett, and Wanner 
1991, Hairer and Wanner 1993). It is also accurate to order 2p + 1 if polynomial of degree p is employed 
to approximate the solution. A key advantage of using the right Radau points (as opposed to the left 
Radau points employed by LaSaint and Raviart 1974) is that the method is simpler and involves fewer 
interpolations in the solution procedure. More importantly, the explicit and implicit schemes here are 
closely related since they use the same intermediate time levels, and the former can serve as a key 
building block in an iterative procedure for the latter. 

The DG derivation here also differs from those in the literature: the standard DG formulation is highly 
algebraic, and the behavior of the solution is hidden in a system of equations; the current formulation is 
geometric, and the behavior of the solution can be observed. 

In addition to the time stepping methods, since the numerical tests have discontinuities, a limiting 
technique that can suppress oscillations near discontinuities while preserving accuracy near extrema is 
discussed. Finally, preliminary numerical results are shown. 

This paper is essentially self-contained and organized as follows. In Section 2, we review Van Leer’s 
scheme III except that our derivation is for arbitrary order. The moment schemes, which have a CFL 
condition of 1 and can be considered as explicit space-time DG methods, are presented in Section 3. 
Extension to the case of two spatial dimensions is discussed in Section 4. Fourier analyses are carried out 
in Section 5 for the one- as well as two-dimensional cases. DG time discretization in the FR/CPR 
framework is presented in Section 6. Section 7 deals with limiting. Preliminary numerical solutions are 
shown in Section 8. Finally, conclusions and discussions can be found in Section 9. 

2. Exact Time Evolution and Projection 

In the fourth installment of a celebrated series of five papers entitled “Towards the ultimate 
conservative difference scheme”, Van Leer (1977) introduced five schemes for advection; among them, 
the first three are piecewise linear and, the last two, piecewise parabolic. Concerning the piecewise linear 
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methods, scheme I is the least accurate but, due to its ease of extension to systems, becomes the most 
popular and is typically known as the MUSCL scheme (monotone upstream-centered schemes for 
conservation laws). As shown in (Huynh 2003), scheme II formulated for a staggered mesh results in the 
CE/SE scheme with e = 1/2 (conservation element/solution element, Chang 1995). Scheme III is 
remarkable since it is third-order accurate in spite of being piecewise linear. The property of having 
higher than expected order of accuracy is called super accuracy or super convergence. The last two 
schemes are piecewise parabolic. Here, scheme IV leads to the standard PPM (piecewise parabolic 
method). Scheme V stores and updates both the cell averages and interface values. Schemes III and V 
have the same error, which is only 1/9 that of scheme IV. As mentioned above, when trying to extend the 
more accurate schemes to the Euler equations, Van Leer “ran into insuperable difficulties”. However, 
scheme III plays a critical role in serving as a guide for the moment scheme (Huynh 2006), which can be 
extended to systems with relative ease. Below, we review Van Leer’s approach for scheme III, which 
amounts to a projection after an exact shift to account for time evolution. Our presentation of the scheme 
is for arbitrary order, not just piecewise linear. 

Consider the scalar advection equation 


u t + au x = 0 


( 2 . 1 ) 


with initial condition 


u(x,0) = u 0 (x) (2.2) 

where t is time, x space, and a the advection speed, a positive constant. By assuming that u 0 is periodic 
or of compact support, boundary conditions are trivial and are omitted. The exact solution at time t is 
obtained by shifting the data curve to the right a distance at, 

“exactO. t) = U 0 (x - at). (2.3) 


Legendre Polynomials. To approximate the solution, we need some preparations. Let Zc be a 
nonnegative integer. Denote by P k the space of polynomials of degree < k. As is standard in the finite- 
volume community, for the reference interval, we use / = [—1/2, 1/2]. Let the inner product of any two 
functions v 1 and v 2 on / be defined by 


On 



Vi(0 v 2 (0 chf, 


and the L 2 norm of a function v by 


INI 



(2.4) 


(2.5) 


Let L k be the Legendre polynomial of degree k defined by applying the Gram-Schmidt process to 
orthogonalize the monomials ... together with the condition that L k (l/2) = 1. The first few 

Legendre polynomials on / = [—1/2, 1/2] are 

L 0 (O = l, Li(0 = 2<f, L 2 (,0 = 6<f 2 — 1/2, 

L 3 (O = 20<; 3 -3$, L 4 (0 = 70<f 4 - 15^ 2 + 3/8. 
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Their plots are shown in Fig. 2.1. Note that compared with the standard Legendre polynomials P k defined 
on [-1, 1], for in [-1/2, 1/2], L k (0 = P k (2ft. 


u 



Fig. 2.1 Legendre Polynomials 


£ 


Discretization. Let the domain of calculation be divided into non-overlapping cells or elements 
denoted by Ej = \xj_ X j 2 , Xj +1 / 2 \ of centers Xj and widths hj. For each Ej, let the local coordinate ^ on 
/ = [-1/2, 1/2] be 

f = (x-xj)/hj. (2.7) 

Next, let p > 0 be an integer. On each cell Ej, we wish to approximate the solution by a polynomial of 
degree p denoted by w ; -. To this end, we need to define the projection operator. Consider the p + 1 
Legendre polynomials L k ,k = 0, ... , p. With t fixed, the solution u(x, t) can be approximated on Ej by its 
projection onto P p : 


v 

Wj (0= £u y , k L k tf) (2.8) 

k = 0 

where the coefficients are given by 

r 1/2 

Uj, k = u{ Xj + fhj, t ) L k ( O dt; / IILJI 2 . (2.9) 

7-1/2 

Thus, on Ej, Wj is the least square fit of degree p for u, and since the projection does not increase the 
norm, we have ||w 7 || < ||u||. The collection of all w ; - as j varies forms a piecewise polynomial function 
denoted by, in the physical coordinate, w(x, t). The function w(x, t) is generally discontinuous at the 
cell interfaces (see Fig. 2.2(a)). 

Note that Uj 0 is the cell average value on Ej, and Uj t represents half of the normalized slope so that 
for p = 1, the value at the right interface is Uj 0 + Uj x . In fact, for all p, the value at the right interface is 
Uj 0 + Uj t + — b Uj p , i.e., each Uj k results in a k- th degree approximation to the right interface value. 
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For simplicity but not necessity, assume that the mesh is uniform: hj = h. Denote the time step by At 
and the CFL number aAt/h by a. Then, since a > 0, the CFL condition takes the form 

0 < ct < 1. (2.10) 

That is, in one time step, the wave advects a distance no more than one cell width. It is this property that 
we wish to preserve when extending the method to one-dimensional systems of equations. 

Let the projection of the initial data be carried out and denoted by, on Ej, 

v 

w°(f,0)= ^u° k L fc (0- 

k = 0 


Exact Time Evolution and Projection. Assume that at time t n , the projection Wy(^, t n ) is known; 
i.e., = Uj k are known for all j and all k < p (the superscript n is understood). We wish to calculate 

all li”/! 1 at time t n+1 (the superscript n + 1 is retained). 

Given the piecewise polynomial data w(x, t n ) = w(x) = (wy(x)}, the exact solution at time t n+1 for 
this data is obtained by shifting the curve w(x) to the right a distance aAt as shown in Fig. 2.2(b): 

v(x) = w(x — aAt). 


That is, on £), with the local coordinate ^ given by (2.7), 


Vj(0 = 


w,— - o + 1), 
W/Of - o'), 


if ^ < -l/2 + o 
if f > -l/2 + o ' 


( 2 . 11 ) 


To evaluate the solution, we split the integral into two parts (again see Fig. 2.2(b)); with k = 0, ... , p, 


n+i _ 
u j, k — 


r —1/2 +(T 

r 1/2 

w y _ 1 (^-a + l)L k (0^ + 

Wjtf - o) L k (0 dt 

y-1/2 

J -l/2+o 



( 2 . 12 ) 


Such a solution can easily be evaluated by a mathematical programming package such as Mathematica. 
For the case p = 0 or piecewise constant data, scheme III reduces to the first-order upwind scheme. 
For the case p = 1 or piecewise linear, the solution is 


u 


n+l 

7.0 


— Uj k + o(Uj_ 1 0 + Uy_ 1(1 — Uj 0 — Uj i) + 0 2 (— Uj- 1?1 + Uj ;1 ), 


u^i 1 = Uj' x + <r(— 3Uj_ 1(0 — 3Uj_ ± i — 3Uj 0 — 3Uj x ) + 

o 2 (/3Uj_i q + 6 Uj—i' \ — 3Uj' o) + — 2Uj—~i' i + 2 Uj' i ) 


(2.13) 

(2.14) 


In matrix form, 
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u 


n+l\ 


7.0 \ 

V U ”W 


a 

— 3a(l — a) 

1 — tr 

3a(l — a) 


W u j-i,o\ 

j V W/-i, i ) 


cr(l - cr) 

— a(3 — 6 ct + 2a 2 ) 

-a(l-a) \( u j,o\ 

(1 — a) (1 — 2a — 2a 2 )/ v u j, i / 


+ 


(2.15) 


« Cell j 



u 



0.5 1.0 1.5 2.0 2.5 3.0 3.5* 


(a) Data 


(b) Solution 


Fig. 2.2 Scheme III (a) Piecewise linear data; (b) Solution obtained by shifting the data a distance aAt 
and calculating the moments of the discontinuous function. 


Note the term a 3 above for the case p = 1. For a general p, the function iv ; (f — a) L p (%) is of degree 
2p in as a result, the term of highest degree in the solution is a 2p+1 . It will be shown by Fourier 
analysis that the method is accurate to order 2p + 1 and is stable for 0 < a < 1. 

Also note that the method is clearly energy-stable: the L 2 -norm of the solution at time t n+1 is bounded 
by that of the data at time t n since the solution is obtained by an orthogonal projection. 

As the final remark of this section, we discuss the extension of this approach. For the advection case, 
the discontinuity at an interface evolves via the exact shift operator. In the case of systems such as the 
Euler equations, in one spatial dimension, these discontinuities result in some combination of a shock, a 
contact, or a fan. Tracking these waves accurately is extremely difficult. Resolving these waves for the 
multi-dimensional cases appears to be an impossible task. 

3. Moment Scheme for the Case of One Spatial Dimension 

Formulation. To avoid the exact shift operator — which makes extensions difficult or impossible — we 
take a different path by using a space-time domain and integration by parts together with the observation 
that the space-time fluxes contain the information provided by the shift operator. Consider the scalar 
conservation law 


u t + f x = 0 


(3.1) 
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where t is time, x spatial coordinate, and / = f(ii(x, t)) the flux. With the wave speed a(ii') = df /du, 
(3.1) implies 


u t + au x = 0. (3.2) 

On the cell E = Ej, let 0 be a test function, i.e., a polynomial of degree p. The conservation law in 
weak form is 


u t 

Jr 


(x t t) 0(x) dx + f x 0 dx = 0. 


* 

Jr 


(3.3) 


We will set <p = L k , the Legrendre basis function, later. Thus, the weak form amounts to projecting onto 
P p . As in the DG formulation, by applying integration by parts to the second term, 

^ J u(x, t ) 0 (x) dx + [f ~ j f 4>xdx =0. (3.4) 

The fluxes / for the boundary term [/ 0]^’^ are chosen to be the upwind ones. At each interface, this 

flux is co mm on for the two adjacent cells and allows the data in these cells to interact. Without this 
common flux, there is no interaction, and the resulting method is inconsistent with physics. 

Integrating in time from t n to t* (t* will be set to the intermediate time levels as well as the final time 
t n+1 later), we obtain 


I u(x, t*) 0(x) dx — I 
J E J E 

-fj 

Jtn Jp 


u(x, t n ) 0(x) dx + J [/ upw 0]*j^* 
/ <p x dx dt = 0. 


dt 


(3.5) 


Denote the surface and volume integral by, respectively, 

SI =-£[/upw Hx'l'l dt (3.6) 

and 


VI = f f f <p x dx dt. 

't n 'e 

The problem reduces to the evaluation of SI and VI, 


„ r. 

u(x,t*) 0(x) dx = u(x, t n ) 0(x) dx + SI + VI. 
Jf Jf. 


(3.7) 


(3.8) 


Here, for 0 = L k and t* = t n+1 , the left hand side above in the local coordinate yields the solution 
hj\\L k \\ 2 xtfX 1 whereas the first term on right hand side, hj\\L k \\ 2 nf k — hj\\L k \\ 2 Uj k . 

Consider now the case of advection. So far, the current approach is equivalent to the exact evolution 
approach of the previous section. Indeed, with the piecewise polynomial data w ; (x) at time t n , the exact 
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space-time solution is w(x, t) = w,(x — at), which has discontinuities along the characteristic lines 
joining the points (Xj_ 1 / 2 , t n ) and (xy_ 1 / 2 + aAt, t n+1 ) as shown in Fig. 3.1. The two formulations, one 
by projecting w(x, t n+1 ), and the other by integrating on the space-time cell Ej X [t n , t n+1 ] yield 
identical solutions, which are exact for the data. At first glance, it appears that we have made the problem 
more difficult since the volume integral takes the form 

VI = f [ aucp x dx dt. (3.9) 

Jt n J e 

In the previous section, we only deal with a discontinuity at one point in space; here, we have to deal with 
a discontinuity along a line in a space-time domain. However, as will be shown, the term <p x comes to the 
rescue in a procedure which ignores the discontinuity. 


Exact Time Evolution 



Fig. 3.1. Space-time cells Ej X [t n , t n+1 ] and the paths of discontinuities (dark line segments) for the 
case of advection 


Time Integration. We proceed by discussing time integration. Due to the term a 2p+1 in the solution 
as noted in the previous section, we need a quadrature with a degree of precision of 2 p in time, i.e., it is 
exact for polynomials of degree 2 p or less (note that we only need a degree of precision of 2 p since the 
vertical edge of £)■ x [t n , t n+1 ] accounts for the additional power in a 2p+1 ). Here, we use to the right 
Radau quadrature with p + 1 evaluation points — consistent with the spatial degree of freedom. Another 
reason for using this quadrature is that the DG method applied to time can be reduced to the Radau IIA 
method (Huynh 2009b), which is a collocation method corresponding to this quadrature. In other words, 
the moment scheme formulated here has the same intermediate time levels and is closely related to the 
implicit Runge-Kutta (RK) scheme resulting from a space-time DG discretization (see §6). 

A brief review of the right Radau quadrature and the associated Butcher matrix is in order (we can also 
employ the left Radau quadrature as in LaSaint and Raviart (1974); the resulting method is, however, 
more complex). Set s = p + 1 where s is the number of stages. Using notations from ODE, on the 
interval [0, 1], let the s right Radau points be denoted by q, l = 1, ... , s where c s = 1. Recall that these 
points result in a quadrature formula on [0,1] that is exact for polynomials of degree 2(s — 1) or less. For 
each l, let the corresponding Lagrange polynomials on [0, 1] be denoted by La ; and defined by the 
condition that La ; (c m ) = 0 for all m ^ l and La ; (c ; ) = 1: 
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T — C, 


(3.10) 


La ; 


m 

c l — c m 


■« = n ^ 

m = 1 1 

The plots of these functions for the case s = 2 and s = 3 ( where s = p + 1) are shown in Fig. 3.2. 


u u 




Fig. 3.2 Lagrange polynomials for the right Radau points: (a) Two points c t = 1/3, c 2 = 1, (b) Three 
points c t = 0.155051, c 2 = 0.644949 , and c 3 = 1. 


Next, suppose the values of a function v are known at these Radau points: v m = v(c m ), m = 1, ... , s. 
We wish to obtain a quadrature formula for J -Ci vir^dT, l = 1, ... , s. To this end, the polynomial of degree 
p = s — 1 interpolating v m , m = 1, ... , s, can be expressed as 


= V r> m La m (r). 

4— im=l 


Integrating from 0 to c h we obtain 


r c l rCl 

v{z)dT=) v m \ La m (r) dr. 
4 _i m= i J Q 


(3.11) 


(3.12) 


Denote the Butcher matrix of dimension s X s by A. The entries a tm of A are given by 

a lm = \ La m (r) dr. (3.13) 

These a im can easily be obtained by using a software package such as Mathematica. The above two 
equations imply, for each l, 


I v(r)dT = / c 
Jo 4_i m= i 


a lm V m - 


(3.14) 


Note that for each row l, setting v = 1, we obtain c ( = £m=i a zm- Also note that for l = 1, ...,p, the 
above quadrature has degree of precision of only p. However, for l = s = p + 1, i.e., for the integral from 
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0 to 1 , due to the definition of the Radau points, the corresponding quadrature has a degree of precision of 
2 p. The column vector C = {cj and the matrix A for the case p = 1 and p = 2 are, respectively 



- 1/12 

1/4 


and 


'0.155051' 

0.196815 

-0.0655354 

0.023771 ' 

0.644949 

0.394424 

0.292073 

-0.0415488 

1 

.0.376403 

0.512486 

0.111111 . 


(3.15) 


Returning to the conservation law, for the space-time domain Ej X [t n , t n+1 ], denote the time level for 
each stage l of the s stages of the above right Radau quadrature by 


t n J = t n + C; At. 


(3.16) 


These are the stages of a Runge-Kutta type scheme where t n,s = t n+1 . See Fig. 3.3(a) and (b). 


Space-Time Solution with No Interaction. Next, we discuss an approximation for the solution on the 
space-time cell £)■ X [t n ,t n+1 ] with no interaction (NI) i.e., no upwinding. This approximation denoted 
by wj v/ (x, t) is a polynomial of degree p, which can be a Taylor series around (xy, t n ). The case of 
p = 1 is trivial: w ; (x) provides u x and (3.2), u t . For a general p, in each cell £), at time t n , the solution 
is known and is a polynomial of degree p in x, 

p 

wj (x)=£u, fe L fe (x). ( 3 • 1 7 ) 

k = 0 

From (3.2), we can calculate all the time as well as mixed space-time derivatives up to degree p, i.e., u t , 
u xt> u tt> ■■■ a t { x j> t n ) via the Cauchy-Kovalevsky procedure (Harten, Osher, Engquist, and 
Chakravarthy 1987). The space-time Taylor series expansion around (xy , t n ) up to degree p then yields 
w^'(x,t) on Ej X [t n , t n+1 ]. 

As an alternative to the above Taylor series expansion, we can use the nodal Lobatto points together 
with a RK time stepping with no interaction. For example, consider the case p = 3. See Fig. 3.3(b). On 
Ej, at time t n , since the cubic w ; (x) is known, we can obtain the values at the four Lobatto points via a 
(4 X 4) matrix multiplication. Next, applying the standard 4-stage RK scheme to (3.2) with no 
interaction, we can evaluate the nodal values at time t n+1 . From the known nodal values for u at t n and 
t n+1 , we can calculate the corresponding u t again with no interaction by (3.2). Then, using a Flermite 
cubic interpolation with the two ends at t n and t n+1 , we can evaluate the nodal values accurate to degree 
p = 3 for all intermediate time levels t n ’ 1 = t n + c L At, l = 1,2, 3. Note that this step is inexpensive and 
the equations for either primitive or conservative variables can be employed with the former more 
economical than the latter. 

An observation concerning the case of advection with a > 0 is in order. We focus, for the moment, on 
only one cell, say, Ej with the polynomial data Wy. For the space-time domain Ej X [t n , t n+1 ], along the 
two interfaces (xy-i/ 2.0 and (xy +1 / 2 ,t), the fluxes are identical to those by method of characteristics 
applied to Wy assuming the domain for Wy(x) is (— oo, oo). Indeed, both the method of characteristics and 
space-time Taylor series produce the exact solution for the polynomial data. However, when interaction 
comes into play, whereas along the right interface, the upwind fluxes reduces to the result by method of 
characteristics for cell j, along the left interface, the values by Wy at x y _ 1 / 2 — a c ( At, namely, 
awy V/ (xy_ 1 y 2 , t n,i ) = Wy(x / _ 1 / 2 — a c ; At), are from outside the cell Ej and will be discarded by the 
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upwind step. In other words, along (Xj_ 1 / 2 , t), the upwind procedure sorts out and employs the values by 
Wj- 1- 

Thus, for general conservation laws, via one of the above two procedures (Cauchy-Kovalevsky or 
explicit RK), for all intermediate time levels t n,i of the right Radau type, at each cell interface x ;+1 / 2 , we 
have two values u L and u R (no interaction). See Fig. 3.3. These values in turn yield the upwind flux 

/upw0'9+l/2> £ ’ )■ 



(a) p = 1 


(b) p = 3 


Fig. 3.3 Space-time nodal points with space via Lobatto and time via right Radau: (a) p = 1, and (b) 
p = 3. The number of spatial points, p + 1, is the same as the number of Legendre basis functions. Along 
each interface, at each Radau intermediate time level, we can obtain the values u L and u R , which in turn 
yield the upwind flux. 


Moment Solution. We can now discuss the moment solution. On the space-time domain £)■ X 
[t n , t n+1 ], with t* = t n ' 1 , switching to the local coordinate ^ on / = [—1/2, 1/2], since dx = hjd%, (3.8) 
yields 


hj\\L k \\ 2 uj;l = hj\\L k \\ 2 Uj k 


£ n , l £ n , l 

~ [ [/upwM-1/2 dt + [ [f(L k )fd(dt. 

J t n J t n J ! 


(3.18) 


For the above last two terms, integrate in (normalized) time from 0 to c L via the right Radau quadrature 
(3.14), we obtain the following surface and volume integrals: with Legendre index k and time index l, 


SIfc — At y a im[/upw(''9+l/2> £ ’ Lk(pCj+ 1 / 2 ) fupw(?Cj—l/2>t ' ) (^ 7 — 1 / 2 )] (3-19) 

<m=l 


and 

K = At J fK’ tn,m XWO d f- 

That is, 


hi\\L h 


ujk = h i 


\ 2 < k 


- siL +vi 


(3.20) 


(3.21) 
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The solution we are seeking corresponds to time level l = s and Legendre indices k = 0, ... , p. 

The surface integral is straightforward. The volume integral is where the difficulty lies. This integral 
must be carried out in a manner that two contrasting requirements hold: (a) the procedure can be extended 
with relative ease to the case of systems but (b) the space-time discontinuity in Fig. 3.1 is somehow 
accounted for. 

What comes to the rescue here is the fact that the degree of (L k )^ is one lower than that of L k ; as a 
result, we can successively calculate uJ'q, , u ; ,! ' 2 , ..., and uj'^ as follows. 

First, for k = 0, the volume integral vanishes since (L 0 )f = 0- Thi s fact and (3.21) yield uJ"q for all l. 

With uJq known, we can calculate uj’l where k = 1 as follows. Since (2.6) implies (Z^)^ = 2, the 
volume integral VI k of (3.20) involves only the projection of / onto 2 L 0 . Thus, we need to approximate / 
up to L 0 but no higher, which is given by /( u n ' 0 ). The solution uj 1 ’/ can therefore be obtained for all l. 

Next, with u"g and u”j ( known, we can calculate uj 1 ^ where k = 2 as follows. Since (2.6) implies 
(L 2 ) f = 12^ = 6L ls the volume integral VI k of (3.20) involves the projection of / onto 6 L 1 . Thus, we 
only need to approximate / up to L 1 . With u”q and uj'l known, we can approximate u by 
u"q L 0 (O + uft Fi(0> an d / by fpQ f 0 (O + fj 1 1 ki(0 (by whatever procedure that provides / to 
the same degree as u). The solution uj 1 '^ can then be obtained for all l. 

Assuming that u"g , ..., uj k _ 1 are known, we now calculate Since (L fe )^ is of degree k — 1, it 
can be expressed as a linear combination of L 0 , . . ., L k _ 1 . In fact, in general, 

Z Ttl \ 1 771 

L 2 ( and (L 2m+2 )^ = 2 y ^ 2 i+i> (3.22) 

i=0 1(=0 


i.e., 


(^ 2 m+i)f — 2L 2m + 2 L 2m -2 + •" and (L 2m+2 )^ — 2L 2m+1 + 2L 2m _ 1 + ••• 

Thus, the volume integral Vlj. by (3.20) involves projecting / onto P k _ ls and we only need to 
approximate / up to L k _ 1 . Since the solution can be approximated by u. n, 0 L 0 (^) + — h 
nfk-i^k- i(f), the flux / can be approximated by i 0 (<f) + — I- //V- i^fc-i(f) ( a g a i n 5 by 
whatever procedure that provides / to the same degree as u) and the solution follows for all l. This 
completes the algorithm. 

Note that in the calculation for uj 1 ' , namely k = p, the last Legendre component k for the above 

successive calculation, we only need to evaluate u the quantities uj’^ where / = 1 , ... , p are not 

needed (we would need them if we use the left Radau or Gauss quadrature or if we use the current 
procedure as a building block in an iterative procedure for the implicit right Radau method). 

Also note that if u is approximated by Wj = Uj 0 L 0 (O + — b u j,k - 1 L k _ i(0, there are several ways 
to approximate / up to L k _ 1 . One way is by evaluating f m = /(wy(f m )) where are the k Gauss 
points and then projecting the polynomial interpolating {/ m } onto P k _ t . 

For the case of advection, the current moment approach is equivalent to the exact evolution approach 
of the previous section. This fact has also been verified by the author via a Mathematica program. 
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For the piecewise linear case, the moment scheme is simple: wj v/ (x, t) is a linear function where, with 
^ on [—1/2, 1/2] and ron [0,1], 


Wj NI tf,T) = Uj 0 + 2 Uj ^ — 2uj t or. 

The extension of the piecewise linear case to the Euler equations is also simple (Huynh 2006). 


Algorithm. The moment scheme is summarized below. Assume that the data Uj k at time t n are given 
for all j and k = 0, ... , p. 

(1) In each space-time cell Ej X [t n , t n+1 ], obtain Wj NI with no interaction by a Cauchy-Kovalevsky 
procedure or an explicit RK with no interaction. 

(2) At each interface Xj_ 1 / 2 , for each right Radau intermediate time level, calculate the left and right 
values using and wj v/ . Obtain upwind fluxes from these two values. 

(3) Successively evaluate u ; n uj 1 / , uj 2 , ..., and u ; n ’p via (3.21) by calculating the surface integral 
via (3.19) and volume integral via (3.20). This completes the algorithm. 

We conclude this section by the following remarks. 

If the upwind fluxes at the interfaces are given (e.g., in an iterative procedure for the implicit right 
Radau time stepping), the above algorithm provides the solution for the space-time cell Ej X [t n ,t n+1 ]. 
This observation shows that the moment scheme can serve as a key building block in an iterative 
procedure for the implicit time-stepping method. 

One way to apply the moment approach to the nodal type of schemes is to transform the nodal values 
into the modal (Legendre) components, carry out the moment calculations, and then transfer back. The 
various transferring, however, involve some costs. The question is whether it is possible to use an 
iteration procedure in the nodal framework with no transferring. 

4 Moment Scheme for the Case of Two Spatial Dimensions 

For the case of two spatial dimensions, the approach of exact time evolution and projection of Section 
2 can be extended in a manner that if the flow is along the diagonal direction, say, from southwest to 
northeast, then the solution recovers the exact data of the southwest cell for the CFL number (a x , a y ) = 
(1,1). Such an approach must allow for interaction among cells that share a common edge as well as cells 
that share a common comer. Its key drawback is that extension to systems is extremely difficult. 
Therefore, the extension of the moment scheme here allows interaction only among immediate neighbors, 
i.e., cells that share a common edge. The price to pay for this simplification is a reduction in CFL 
condition when the flow is not along the edge direction. 

To proceed, with u = u(x, y, t ), the 2D conservation law takes the form 

t if fx "L 9y 0- 

Let initial condition be u(x, y, 0) = Uj nit (x,y) . For the Euler equations, u , /, and g are vectors of four 
entries. With no loss of generality, we consider only the scalar case. Assume that the solution u is 
periodic or of compact support so that boundary conditions are trivial. 
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Denote / = (/, g ) and V = Then V ■ / = f x + g y , and the above can be written as 

u t +V-f= 0 . ( 4 . 1 ) 

Let the domain of calculation be divided into a non-overlapping unstructured mesh of cells Ej which 
are quadrilaterals or triangles. On each cell Ej, let the solution u be approximated by Uj of degree p. If the 
cell is a quadrilateral, we can use either a polynomial in ^ and g of degree p or less, or we can use tensor 
products of degree p x p. For example, if p = 1, we can use either a linear polynomial a + b<f + eg, or 
for the case of tensor products, a bilinear polynomial a + b<f + eg + d^g. 

Using Uj, the fluxes / and g can be approximated by fj and g j of the same degree and type as Uj. The 
piecewise polynomial functions { Uj }, {/)■}, and {g j} can be and usually are discontinuous across cell 
interfaces. 

From here on, unless otherwise stated, we focus only on the approximate polynomials. For simplicity 
of notation, the subscript j is understood: Uj, fj , and Ej are abbreviated to u, f, and E when there is no 
ambiguity. 

Let 0 be a test function on E, i.e., a polynomial of the same type as u. We require the solution to 
satisfy, for the moment, formally, 



dE + 



dE = 0 . 


As in the ID case, integrating by parts and applying the upwind flux for the boundary term to account for 
data interaction, the solution is required to satisfy 



L 


dE + I [f ■ 0 ds 



dE = 0 . 


( 4 . 2 ) 


Integrating in time from t n to f * = t n ’ 1 , 



,y, t*) ip dE 



,y, t n ) <p dE — SI - VI 


= 0 


( 4 . 3 ) 


where the surface and volume integral are, respectively, 


and 



C f ^)upw 0 ds dt 


VI 



Vcp-fdE dt. 


( 4 . 4 ) 


( 4 . 5 ) 


The problem reduces to the evaluation of SI and VI. 

To this end, as in the ID case, we first obtain w , j v/ with no interaction by a Cauchy-Kovalevsky 
procedure or an explicit RK with no interaction up to appropriate order. For the case of tensor products, 
we would need to derive time derivatives of degree up to 2 p if the Cauchy-Kovalevsky procedure is 
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employed. Next, along each edge, for each right Radau intermediate time level, calculate the left and right 
values using appropriate Wj" . Obtain normal upwind fluxes from these two values. 

The surface integral again poses no difficulty. As for the volume integral, we make use the fact that 
both components of V0, namely, (p x and (p y are of at least one degree lower than (p. As in the ID case, 

we can evaluate the volume integral and successively calculate the Legendre components u, n ' 0 l 0 , uj’l 0 , 
u"p i, uj ‘2 0 , u?l lt uj 02 , ... This completes the algorithm. 

5 Fourier Analyses 

One-Dimensional Case. For the ID case, consider the advection equation with positive speed a, 

u t + au x = 0. 


The cells are Ej = [/' — 1/2, j + 1/2]. Denote by Uj the column vector of K = p + 1 components, 

T 

Uj (llj 0 , tty 1, ... , Uj p) 

where T represents the transpose. Let At be the time step and set a = aAt. Assume that the data U 1 - = Uj 
are known. The solution Uj +1 can be expressed as 

U? +1 = C 0 Uj + C_ 1 Uj_ 1 (5.1) 


where C 0 and are K X K matrices. As examples, for p = 1, i.e., the linear case 

)• c » = (- 


a er(l — a) 

3ct(1 — a) —a( 3 — 6 o + 2a 2 )/ 


C-i= ( 

and, for p = 2 or the parabolic case, 

C - 1 = 


l — o cr( 1 — a) > 

3a(l — a) (1 — <r)(— 1 + 2a + 2a 2 )/ 


a 


a(l — a) a( 1 — a)( 1 — 2a) 

—3a(l — a) —a(3 — 6a + 2a 2 ) —3a(l — a)(l — 3a + a 2 ) 

5er(l — cr) (1 — 2cr) 5cr(l — cr) (1 — 3 ct + a 2 ) a( 5 — 30ct + 50cr 2 — 30cr 3 + 6CT 4 ); 


and 


C 0 = 


l — o —o(l — o) —a(l — a)(l — 2a) 

3a(l — a) (1 — a)( 1 — 2a — 2a 2 ) —3a{l — cr) (1 — a — a 2 ) 

K — 5ct(1 — a)( 1 — 2a) 5ct(1 — cr) (1 — a — a 2 ) (1 — cr) (1 — 4cr — 4 ct 2 + 6 a 3 + 6cr 4 ), 

Note the highest degree term for p = 1 is a 3 and, that for p = 2, cr 5 . 


Stability. Let the wave number be denoted by w, - n < w < 7r, and the imaginary number, by / 
instead of the standard notation i to avoid confusion with the cell index i later. Assuming the solution is a 
harmonic such that Uj = e~ Iw Uj , we obtain the amplification matrix 

A = C 0 + e~ IW C_ 1 (5.2) 


which satisfies, by (5.1), 


NASA/TM— 2013-218077 


16 



(5.3) 


Uf +1 = AUj. 

With each value of a and w, the corresponding matrix A has p + 1 = K complex eigenvalues. For 
stability, these eigenvalues must have magnitude no larger than 1. Numerical calculations for these 
eigenvalues for all p tested by this author result in stable schemes provided that 0 < a < 1. 

Figures 5.1 show the magnitude of the amplification factors (or eigenvalues) for the case p = 1 and 
Figs 5.2, the case p = 2. 




Figs. 5.1 Magnitude of the two amplification factors (or eigenvalues) for the case p = 1. 



(a) Principal 


(b) Spurious 


(c) Spurious 


Figs. 5.2 Magnitude of the three amplification factors for the case p = 2. 


Accuracy. To calculate the order of accuracy, denote the principal eigenvalue of A by e p = e p (a,w ). 
This value approximates e~ Iaw . A scheme is accurate to order m if, for a fixed a, e p approximates e~ I<JW 
to 0(w m+1 ) for small w, 


e p — e~ Iaw = 0(w m+1 ). 


(5.4) 
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In practice, it is difficult if not impossible to derive a Taylor series expression for e p when p > 1. 
Therefore, we obtain the order of accuracy of a scheme by a numerical calculation. First, set a = a 0 = 
0.8, say, and w = w c = 7r/4. We can calculate the coarse mesh error 

er c = e p Oo , w c ) - e _,£r ° Wc . (5.5) 

By halving the wave number, Wf = w c /2 = 7r/8 (equivalent to doubling the number of mesh points), the 
corresponding fine mesh error is 


ei/ = e p {a 0 ,w f ) - e Ia ° w f. 
For a scheme to be m - th order accurate, after one time step, 


er c 

er f 


2 m+1 


That is, 


m 



Log(2) 


- 1 . 


(5.6) 


(5.7) 


(5.8) 


Thus, for a scheme to be of order m, when we march to a certain final time, doubling the mesh results in 
an error smaller by a factor of 2 m . 

It can be verified, as will be shown in the more involved 2D cases, that a moment scheme using 
polynomials of degree p is accurate to order 2 p + 1. 


Two-Dimensional Case. We next briefly describe Fourier (Von Neumann) analysis for the 2D case. 
On the domain (-00,00) x (-00,00) , consider the advection equation 


u t + au x + bu y = 0. 


The cells are the squares = [i — 1/2, i + 1/2] x [/' — 1/2, y + 1/2] centered at (1,7). For the case of 
polynomials of degree p, the number of Legendre basis functions is K p = (p + l)(p + 2)/2. For the case 
of tensor products, the number of Legendre basis functions again denoted by K p , is K p = (p + l) 2 . In 
addition, denote by l/y the column vector of K p components, 


Uij — (ttjj.o.oi u i,j,i,o> u i,j, o,i’ ■■■ ) 


(5.9) 


As an example, for p = 1, 


j 

Uij = ( u i,j, 0,0- u i,j, 1 , 0 ’ u i,j, 0 , 1 ) 


Or, for the case of tensor products, 


T 

Uij = ( u i,j, 0,0- u i,j,l,0> u i,j, 0,1 1 u i,j,l,l) ■ 
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Set a x = aAt and a y = bAt. The 2D moment solution can be written as 


Uij 1 = Co.oUij + C_ lj0 y t _ w + 


(5.10) 


Here, C 0 0 , C_ 10 , and C 0 ,_ 1 are K p X K p matrices whose entries involve er x and a y . As an example, for 
the linear case (5.9), 


G Y 


3(J x (Jy 


~O x Oy 


C o — g x ) 


6(7^. H“ 2(J % ) 2(7^.) ] t 

CT X (7y(3 2(7^.) 


°x(l - 2cry) 


and 


C 0 .-i — 


(Tv 


~G X Oy 


3(J X (Jy 


Gy(l ~ 20£) 

~ 3 Gy ^ 1 (Jy ^ Oy ^ 3 2 Oy ^ 


(Ty ( 1 (Ty ^ 

G X Gy(3 2 (Ty^ 

-cr y (3 - 6er y + 2cr y ) / 


(T v 


G-\j 


Q,0 — 


- 3 cr :r (— 1 + er* + cr y ) 1 — 3 er x + 2 er| — a y + 2 a x a y 

\ 3(Ty( 1 + G X + (Ty) 


®jjc( 1 3" 0 X + (Ty) 

+ : 

20jj[ .Gy^Px b @y) 


(Ty ( 1 + (T X + (Ty) 

20 x Oy{a x "F Oy) 
Ty + 2(T x (Ty 


1 — (T r — 3 (T v + 2 (T r (Tv + 2 (Tv 


Stability. Let the wave number in the x -direction be denoted by w x and in the y-direction, iv y 
(involving no partial derivative), and let the imaginary number be denoted by I. Similar to the ID case, 
replacing Ui-ij and Uij-i respectively by e~ IWx Uij and e~ lw yUij, we obtain the amplification matrix 

A = C 00 + e~' Wx C_ 10 + e~ Iw yC 0 _ 1 . (5.11) 

Consider a fixed value ( a x , a y ). For each wave number ( w x , w y ) where - n < w x < n and - n < w y < n, 
the corresponding matrix A has K p eigenvalues. For stability, these eigenvalues must have magnitude no 
larger than 1 for all wave numbers. The domain of all values ( a x , a y ) satisfying this stability condition is 
called the stability region of the scheme. 

Figure 5.3 shows the stability region for a square mesh and polynomials of degree p = 1, ...,4. As is 
well-known, the CFL limits along the diagonal for the first-order upwind scheme (p = 0) is V2/2 ~ 
0.7071. In the case of the moment schemes, for p = 1, ... ,4, they are, respectively, 

0.4714, 0.2859, 0.2190, and 0.1654. 
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Figure 5.3. Stability regions for schemes using polynomials of degree p. Note that for p = 3 and 4, 
the stability conditions along the x and y directions are 1, but these conditions reduce very quickly when 
the flow slightly deviates from the axes directions. 


Figure 5.4 shows the stability region for a square mesh and polynomials via tensor products of degree 
p X p where p = 1, 2, 3, and 4. The CFL limits along the diagonal for p = 1, ... ,4 are respectively 

0.4714, 0.2834, 0.2526, and 0.1702. 
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Figure 5.4. Stability regions for schemes using polynomials via tensor products of degree p X p. 
Again, for p = 3 and 4, the stability conditions along the x and y directions are 1, but these conditions 
reduce very quickly when the flow slightly deviates from the axes directions. 


Thus, the reduction in the CFL condition for flows along the diagonal direction is significant. The 
main reason is that information does not directly communicate among cells that share only one comer. 
That is, for a flow from, say, southwest to northeast, information flows to the cells to the east and the 
north before it can reach the cell at the northeast corner. This difficulty concerning information flow can 
be remedied by using the moment scheme as a procedure to update the solution in an iterative process for 
an implicit space-time DG scheme as will be discussed in the next section. 
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Accuracy. To calculate the order of accuracy, using the notation of the ID case, denote the principal 
eigenvalue of the amplification matrix A defined in (5.11) by e p = e p (a x , a y , w x , iv y ). This principal 
eigenvalue approximates e -/ (°* w x+tfyw y ) -p 0 fj nt | the order of accuracy of the scheme, first, fix ( o x ,o y ), 
say, ( a x , Oy) = (0.2, 0.1). For the coarse mesh, set (w % , w y ) = (n/4, n/6). We can calculate the coarse 
mesh error 


er c = e p (a x , a y , w x , w y ) - e Ka x w x +a y w y ) (5.i2) 

Next, halving the wave number, set (w % , vv y ) = (7 t/8, n/12). A formula similar to the above yields the 
corresponding fine mesh error er >. For a scheme to be m-th order accurate, (after one time step) 


Log 



Log(2) 


1 


(5.8) 


As can be seen in Tables 5.1 and 5.2, the 2D moment schemes using polynomials either of degree p or 
via tensor products of degree p x p are accurate to order 2p + 1. 


Polynomial Degree 

Coarse Mesh Error 

Fine Mesh Error 

Order of Accuracy 

1 

1.88 X 10“ 3 

1.28 X 10“ 4 

2.87 

2 

2.58 x 10“ 5 

4.26 X 10“ 7 

4.92 

3 

1.83 x 10“ 7 

7.47 X 10- 10 

6.94 

4 

7.22 X 10 -10 

7.49 X 10“ 13 

8.91 


Table 5.1. Errors and order of accuracy of 2D moment schemes using polynomials of degree p; 
(a x , Oy) = ( 0 . 2 , 0 . 1 ), coarse mesh (w % ,w y ) = ( n/A , n/6), and fine mesh (w % ,w y ) = ( 7 t / 8 , n/12). 


Polynomial Degree 

Coarse Mesh Error 

Fine Mesh Error 

Order of Accuracy 

lxl 

8.3 X 1(T 4 

5.42 X 10“ 5 

2.94 

2x2 

3.09 X 10“ 6 

4.83 X 10“ 8 

5. 

3x3 

9.47 X 10“ 9 

3.81 X 10- 11 

6.96 

4x4 

1.14 X IQ -11 

1.41 x 10“ 14 

8.66 


Table 5.2. Errors and order of accuracy of 2D moment schemes using tensor products (degree p X p); 
(a x , er y ) = (0.2, 0.1), coarse mesh (w x , w y ) = (n/4, n/6), and fine mesh (w x ,w y ) = (7t/8, n/12). 
Due to the larger degree of freedom, errors for schemes via tensor products are smaller than the 
corresponding errors for schemes via polynomials of degree p. 
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6 Implicit Space-Time DG Schemes 


The moment scheme discussed in Section 3 is closely related to the implicit space-time DG method of 
this section. When applying the DG approach to discretize time, as shown by LaSaint and Raviart (1974) 
who employed a quadrature of left Radau type, the result is an implicit Runge-Kutta scheme. It was later 
shown in (Huynh 2009b) that using the FR framework, the right Radau points are a natural choice for 
quadrature points. The resulting DG method is simpler and is identical to a collocation method called 
Radau IIA (Lambert 1991, Hairer, Norsett, and Wanner 1991, Hairer and Wanner 1993). We review the 
DG method and apply the FR approach to time marching to show the close relation between the resulting 
implicit space-time method and the moment scheme. 

Consider the ODE 


u'(t) = f(u,t), (6.1a) 

with initial condition 

u(t°) = u°. (6.1b) 

Here, some of the notations are from ODE; as such, the above / is not related to the flux / of Section 4. 
Next, let At be the time step size and set t n = nAt, n = 0, 1, 2, ...; the assumption of constant step size is 
only for convenience; in practice, time step sizes can vary. When dealing with time, we use / = [0, 1] and 
t = t n + zAt where r varies on 7; thus, t varies on I n = [t n , t n+1 \. 

Assume that the data u n at time t n is known, we wish to calculate the solution u n+1 at t n+1 . For 
n = 0, u n is the initial condition u° in (6.1b). 

The DG method seeks to approximate the solution for (6.1) by a function denoted by u h which can be 
and usually is discontinuous across each t n . On 7 n , this function is a polynomial of degree s — 1, which 
is defined by its values at s points at locations corresponding to the stages. Note that when applying these 
time-stepping methods to conservation laws, if p is the degree of the spatial polynomial, the number of 
stages in time is 


s = K = p + 1. 

Due to the property of super accuracy in time discussed later, we could combine parabolic approximation 
in space or p = 2 with linear approximation in time or s = 2. Such combinations remain to be explored. 

Assume that the calculation for 7 n_1 = [t n_1 , t n ] is finished, and we have the solution u h on that 
interval. Thus, the value just to the left of t n is known, namely u h Recall that the solution can be 

and usually is discontinuous across each t n . We assume time marches forward, thus, the value at t n is 
that from the left (small square of Fig. 6.1), 

u n = u h (t n ~). (6.2) 

We wish to calculate u h on I n = [ t n , t n+1 ], i.e., to obtain 

u h (t n+1 -) = u n+1 . (6.3) 


NASA/TM— 2013-218077 


23 



u 



Fig. 6.1. DG solution u h can be and usually is discontinuous across each t n . The value at t n is chosen 
to be that from the left (small squares): u n = u h (t n “). 


DG formulation. On 7" = [t n ,t n+1 ], let 0 fe be a set of basis functions for Pk-i, k = 1, The 
solution is approximated by 


u h = 


K 

^ ' u k 0k- 
k=l 


Let 0 be a test function, which is in P K _ x ; typically, 0 is one of the basis functions 0 fe , k = 1, 

Again on 7 n , the weak form for u' = f is, formally (the correct equation is (6.6) or (6.7)), 

(u',0) = (/,0). 

Let f h be the projection of / onto P K - Since 0 varies on P K -i, we have (/, 0) = ( f h , 0), i.e., we can 
project first, then apply the inner product. The above implies 

(u',0) = (/ h ,0). 


We wish to find u h in P K -x such that for all 0 in Pk-i, again formally, 

«,0) = (/ h ,0). (6.4) 

For the DG method, u h (t n+ ) is allowed to be and is generally different from u h (t n_ ) =u n . By 
applying integration by parts to the left hand side above, 

K^][n +1 - 0/1-0') = (/k-0)- 


Concerning the boundary term, there is no ambiguity with test function 0, but for u h , since it is 
discontinuous across t n and t n+1 , a more precise expression is 

[u h 0]^n +1 = u ft (t n+1 “)0(t n+1 ) -u ft (t n+ )0(t n ). 

To involve the initial condition, however, we replace u h (t n+ ') in the boundary term by u h (t n “) = u n in 
(6.2). Thus, the DG solution u h is required to satisfy, for all 0 in Pk-i, 

u /l (t n+1 )0(t n+1 ) - u h(t n - )0(0 n ) - (u-tv 0') = ( fh. 0)- (6-6) 
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Applying integration by parts again to (u h , 0'), 


-\u h (t n ) - u h (t n+ )] 0(t n ) + «,0) = (6.7) 

Note that this time we use u h (t n + ) for the boundary term u h (t n ). 

Also note that the above, the result of integration by parts twice, is the ‘strong form’ whereas (6.6) is 
the ‘weak form’. 


Flux Reconstruction formulation. Our goal is to eliminate the test function in (6.7) so that the 
integral formulation results in a differential one. To this end, we raise the following question: can we find 
a polynomial ,g LB = g on I n which possesses the property that, for any 0 of degree K — 1 or less, 

— 0(t n ) = Cg',0). (6.8) 

If such a polynomial exists, we can combine ( u' h , 0) in (6.7) with ( g' , 0 ) and eliminate the test function. 
From (6.7), since f h is of degree K — 1, we require g' to also be of degree K — 1; as a result, g is of 
degree K. Applying integration by parts to the right hand side of (6.8), we have 

-0(t") = <7(t n+1 )0(t n+1 ) - 5 (t n )0(t n ) - ( g , 0'). (6.9) 

The above holds if g satisfies, 

g(t n ) = l, g(t n+1 ) = o, (6.10) 


and, for all 0 in P K _ l5 


(< 7 , 0 ') = 0 . ( 6 . 11 ) 

Since 0 is of degree K — 1, 0' is of degree K — 2; moreover, 0' spans P K - 2 as 0 spans P K _ Eq. (6.1 1) 
implies that g is orthogonal to P K _ 2 , i.e., for any polynomial cp of degree K — 2, 

(g,cp) = 0. ( 6 . 12 ) 

Orthogonality to P K _ 2 provides K — 1 conditions; (6.10) provides the other two. 

Let the right Radau polynomial of degree K be denoted by R r k and defined by the conditions that 
Rr,k (t n ) = 1 and R r k vanishes at the K right Radau points on I n . See Fig. 6.2(a). We wish to show that 

9 = Rr,k- (6.13) 

That is, (6.10) and (6.12) hold with g replaced by R r k . Indeed, (6.10) is immediate by the definition of 
Rr.k- Next, since R r k is of degree K , for any <p of degree K — 2, R Ri k<P is of degree 2(K — 1), and 
R r , K qj vanishes at the K right Radau points. In addition, the right Radau quadrature is exact for R r K <p. 
Consequently, (6.12) holds, and (6.13) follows. 

The answer to the question posed for (6.8) is thus positive. Note that g = ,g LB , which corrects for the 
jump at the left boundary, is the right Radau polynomial (vanishing at t n+1 ). Loosely put, the condition 
g(t n ) = 1 deals with the jump at the left boundary whereas (6.12) together with g(t n+1 ) = 0 implies 
that away from the left boundary, g approximates the zero function. 
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We now return to our task of eliminating 0. With g defined above, (g' , 0) = — 0(t n ). Consequently, 
by (6.7), 

[u h (t n ~)-u h (t n+ )] (g 1 , 0) + (Uft,0) = (/ h ,0). (6.14) 

What is crucial here is that 0 can be factored out and canceled and, since u n = u h (t n “), 

([u n -u h (t n+ )]g + u h y = f h . (6.15) 

Therefore, [ u n — u h (t n + )]g + u h is a polynomial of degree K whose value at t n = t n+ is, by (6.10), 

[u h (t n ~) -u h (t n+ )]g(t n ) + u h (t n+ ) = u h (t n ~ ) = u n . (6.16) 


Thus, to deal with the jump [u h (t n ) — u h (t n + )] = u n — u h (t n + ) at t n , we set 

U = [u h (t n ~) -u h (t n+ )]g + u h . (6.17) 

See Fig. 6.2(b). Differentiating U and employing (6.15), we obtain 

U\t) = f h {t,u h (t)). (6.18) 

The function U is the key idea of FR. On [t n , t n+1 ], assume that u h which does not match the 
boundary condition u n is given. We can add the correction term [u 11 — u h (t n+ ')]g to u h as in (6.17) so 
that the resulting U takes on the left boundary value u n , does not alter the value u h (t n+1 ~) = u n+1 , and 
retains the property of ‘best possible’ approximation to u h . The function U is continuous across cell 
boundaries (as opposed to u h which is discontinuous). To put it differently, u h is of degree p; U is of 
degree s = p + 1 determined by p + 2 conditions; two conditions are known, namely, U (t n ) = u n and 
t/(t n+1 ) = u h (t n+1 “); the remaining p = s — 1 conditions are given by the requirement that the 
projection of U and u h onto P v -\ are the same. This last requirement is the result of the fact that g is 
orthogonal to P p -\ by (6.12). 

Denoting the K right Radau points on I n by t n,k , k = 1, ...,K where t n,K = t n,s = t n+1 . See also Fig. 
3.3. Since g vanishes at these points, (6.17) implies 

U(t n ' fe ) = u h (t n ’ k ). (6.19) 


By employing these Radau points as collocation points (or quadrature points), (6. 1 8) leads to 

U'{t n ’ k ) = f(t n ' k , U n ’ k ). (6.20) 

Here, we approximate f h , which is of degree p = K — 1, by using the K values f(t n,k , U n,k ), k = 
1 K. 

In summary, to calculate the DG solution, we first obtain the polynomial U by solving the implicit 
equations (6.20); here, U is of degree K = s interpolating the K + 1 values at t n and at the K right Radau 
points: t/(t n ) = u n and I/(t n,fe ) = U n,k , k = 1, The DG solution u h of degree K — 1 is defined by 
the values U at these K right Radau points. 

The scheme defined by (6.20) is identical to an implicit collocation Runge-Kutta scheme using the 
right Radau points and accurate to order 2p + 1 = 2 K — 1 called Radau IIA (Lambert 1991, Hairer, 
Norsett, and Wanner 1991). 
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0.0 0.2 0.4 0.6 0.8 1.0 


(a) Radau polynomial R r 3 



Figure 6.2. (a) The correction function g, i.e., the right Radau polynomial of degree s = 3 on 
/ = [0,1]; (b) u h is of degree s — 1 = 2 and discontinuous across t n ; U is of degree s = 3, continuous 
across t n , and defined by: U(t n ) = u n , U(t n+1 ) = u h (t n+1 and the projections of U and u h onto 
P s —2 (providing s — 1 conditions) are the same. 


We now discuss the relation between the current implicit time stepping and the moment scheme 
(explicit time stepping) of Section 3. For the conservation law u t + f x = 0, with the same spatial 
discretization using the Legendre polynomial, the key difference between the two time stepping methods 
is the following. For the moment scheme, the Cauchy-Kovalevsky procedure, i.e., the space-time Taylor 
series expansion with no interaction denoted by w. NI applies. For the implicit scheme, this expansion must 
be replaced by the solution uj 1 '^ itself; recall that the Legendre components correspond to k = 0, ...,p, 
and the time stages, l = 1, ... , s. As such, the moment method of updating the solution can be employed in 
an iterative manner to calculate the implicit solution. This subject, however, is beyond the scope of the 
current paper. 


7. A Limiting Technique 

As is well-known, high-order methods generate oscillations near a shock or a discontinuity. One way 
to avoid oscillations is to limit the data so that it is not ‘too steep’ before updating. The problem with 
limiting is that it can cause a loss of accuracy near extrema. We explore ideas that that can suppress 
oscillations near shocks while preserving accuracy near extrema here. The limiting procedure below is 
designed for the piecewise linear case with one spatial dimension. The procedure employs the techniques 
introduced in (Huynh 1995) and (Suresh and Huynh 1997) but has a stencil of only three cells and is 
simplified by requiring that the magnitude of the slope lies in a certain range. 

To proceed, assume that Uj 0 and u< i are given for all j. We abbreviate the cell average values: 

tt/,0 t ij. 

As discussed after (2.9), the first moment u,- 1 represents half of the normalized slope so that the value at 
the right interface of cell j is Uj 0 + Uj :1 . Next, set 

X/ Uj Uj _ 1; s R Uj+ 1 Uj, (7.1) 


and 
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s L imi = 2min(|s L |,|s R |). 


(7.2) 


The limit s Liml for the normalized slope is similar to Van Leer's limit (1977) except that we use the 
magnitudes. It can suppress oscillations near a discontinuity, but has the drawback of causing a loss of 
accuracy near extrema. With 2 1 Uj 1 1 representing the normalized slope, if 

2 Ki| ^ s Liml, (7.3) 

limiting does not alter the slope, and we move on to the next j index. At most of the smooth regions of the 
flow, the above condition is satisfied; thus, limiting is typically needed only for a small number of cells. 

If (7.3) does not hold, the cell is near a discontinuity or an extremum, and the following calculations 
are carried out. Our goal is to enlarge the limit in a manner that near a shock, it reduces to (7.1), but near 
an extrema, it is larger and provides enough ‘room’ so that the original slope is not altered. To put it 
differently, at smooth regions, the results by the moment or DG schemes are highly accurate and should 
not be altered. 

As will be shown by numerical examples, the moment scheme produces, for each cell, a slope that is 
downwind biased if a is near 0, and upwind biased is a is near 1. Therefore, we take extra care to account 
for this slope variation. 

For any real numbers x and y, let minmod(x, y) be the median of x, y, and 0, 

1 

minmod(x,y) = -[sign(x) + sign(y)] min(|x|, |y|). 


At interface j — 1/2, set 

Vu = ( uj + Uj_ 1 )/2, v 12 = uj_ t + Uj_ x 1 - v lv v 13 = uj - u jj j-%. (7.4) 

Loosely put, we use v tl as a pivot; v 12 is from the left, and v 13 , the right; typically, near a smooth 
extrema, v 12 and v 13 are of the same sign and relatively close to each other. Set 

v 14 = minmod(u 12 , v 13 ). (7.5) 


At interface j + 1/2, set 

V 21 = ( U ; + U j + 1)/2< ^22 = U j + 1 — U j+ 1,1 — V 21< V 23 = U ] + U j, 1 — ^215 


and 


v 24 = minmod(u 22 , ^> 3 )- (7-6) 

Near a smooth extremum, v 14 and v 24 are typically of the same sign and relatively close to each other. 
For the cell j, set 


Vmm = |minmod(u 14 , u 24 )|, v mx = max(|u 14 |, |u 24 |). (7.7) 

Thus, v mm and v mx are also relatively close to each other near a smooth extremum. Near a discontinuity, 
however, they are far apart. With ‘ra’ short for ratio and ‘ap’ for accuracy preserving, set 
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v ra ~ , -t n— 20 ’ Va P ~ 6v mm V ra . (7.7) 

I/ mx ' xu 

Note that v ra is a non-dimensional quantity. Near a discontinuity, it is close to 0; near an extremum, it is 
relatively close to 1. Therefore, v ap provides ‘room’ for the interface value near a smooth extremum. 
Next, set 


s min l< l s R I) < s Lim2 s min 2v ap , s Lim3 max ( s Liml< ^Limz)- (7-8) 

Note the factor 2 accounts for transferring a quantity of type v to type s. Finally, limit Uj x by 

uj, i «- sign(u y<1 ) min(|u ;jl |, s Lim3 /2). (7.9) 

Note that if (7.3) holds and we still apply the above limiting procedure, the result recovers the original 
Uj x . Thus, (7.3) only serves the purpose of saving computing time; whether it is applied or not, the final 
solution is the same. 


8. Numerical Results 

The following preliminary numerical tests show the results by the piecewise linear moment (or upwind 
moment) scheme. Recall that for ID advection, the scheme is identical to Van Leer’s scheme 3. 

The first test is for an advection with constant speed 1. The initial condition is a sine wave on the 
domain [0, 1] with 6 cells and periodic boundary conditions. The continuous (red) curve represents the 
exact solution, the dots the cell average numerical solutions, and the line segments, the piecewise linear 
solutions. Figs. 8.1 and 8.2 show the following behavior of the scheme: for each CFL number, the scheme 
has a ‘preferred mode’ or eigenvector. This eigenvector corresponds to a piecewise linear function with 
slopes that are downwind biased if a is close to 0, and upwind biased if o is close to 1 . 

Figs. 8.1 shows the results with a = 0.1. For small er, the solution Tikes’ downwind biased slopes in 
the sense that no matter what kind of slopes we start with, the solution will turn them into downwind 
biased slopes after a few time steps. In (a), for each cell j, we start with slopes biased in the ‘wrong’ 
direction, i.e., upwind biased slopes defined by the exact interface value Uj _ t / 2 and exact cell average 
value Uj. Fig. 8.1(c) shows the solution after 10 time steps; the wave travels a distance of one cell. At the 
maximum, the slope is now slanted downward, opposite to the initial slope in (a). In fact, in (b), after just 
4 time steps, the solution already exhibits this ‘preferred mode’ for the slopes: downwind biased for small 
a. Fig. 8.1(d) shows the solution after 60 time steps; the wave advects 1 period. Comparing (a) and (d), 
again note the switch in bias of the slopes. 

Figs. 8.2 shows the results with a CFL number of 0.9677, which is the exact opposite compared with 
Figs. 8.1. For CFL number close to 1, the solution Tikes’ upwind biased slopes. We again start with 
slopes biased in the ‘wrong’ direction, i.e., downwind biased slopes defined by the exact cell average 
value Uj and exact interface value Uj +1 / 2 as shown in Fig. 8.2(a). After 31 time steps, the wave travels a 
distance of 30 cells, i.e., the profile advects 5 periods. The solution is shown in Fig. 8.2(b), where the 
slopes are upwind biased. 

The above property of biasing toward the upwind or downwind sides holds true for higher-order as 
well (p > 1). In that case, the p + 1 right Radau and left Radau points come into play in a manner similar 
to the linear case above. (The upwind biased linear functions in Fig.8. 1(a) interpolate the values at the two 
left Radau points for each cell, and the downwind biased ones in Fig.8. 2(a), the values at the two right 
Radau points.) 
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Fig. 8.1. CFL number a = 0.1. The continuous (red) curve represents the exact solution, (a) Piecewise 
linear data with slopes biased in the ‘wrong’ direction, i.e., upwind biased (slanted upward near the 
maximum); (b) after just 4 time steps, solution exhibits a ‘preferred mode’ for the slopes, which is 
downwind biased since a is close to 0; (c) solution after 10 time steps; the wave travels a distance of one 
cell; and (d) solution after 60 time steps; the wave advects 1 period. 


y y 




Fig. 8.2. CFL number a = 0.9677, number of time steps 31, distance of advection 5 periods. The 
continuous (red) curve represents the exact solution, (a) Piecewise linear data with slopes biased in the 
‘wrong’ direction, i.e., downwind biased, or slanted downward near the maximum; (b) the solution 
exhibits a ‘preferred mode’ for the slopes, which is upwind biased for a close to 1. 
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Note that the solution in Fig. 8.1(d) is damped more than that in Fig. 8.2(b) in spite of the fact that the 
wave travels much further in the latter. The reason is that for Fig. 8.2(b), the error associated with the 
CFL number of 0.9677 is small (er = 1 has no error), and the number of time steps is only half that of 
Fig. 8.1(d). In other words, the errors accumulate more in Fig. 8.1(d) due to the larger number of time 
steps. Since a = 1 has no error, for the moment scheme (as opposed to the case of RK time stepping), to 
minimize error, we take time steps as large as possible. 

Note that if we use piecewise linear DG with RK time stepping, the solution behaves in a manner 
similar to that of Fig. 8.1, i.e., it Tikes’ downwind biased slopes. The moment scheme, however, has this 
peculiar property of shifting the slopes it Tikes’ to upwind biased when a is close to 1. It appears difficult 
to devise a semi-discrete formulation with this shifting ability. 

The second test is again advection. The initial condition consists of a Gaussian, a rectangular, and a 
semi-ellipse wave on the domain [0,1] with 100 cells and periodic boundary conditions. The CFL 
number is .8, and the final time is t = 1, i.e., the profile advects 1 period after 125 time steps. The 
solutions by the standard upwind scheme using the Van Albada limiter and the upwind moment scheme 
using the limiting procedure in Section 7 are shown in Fig. 8.3. The dots represent the cell average values. 
Note that the standard upwind solution in (a) has a loss of accuracy (only first-order accurate) near 
extrema, and discontinuities are smeared. The upwind moment solution has no loss of accuracy near 
extrema, and discontinuities are resolved with four points. 




(a) Standard upwind scheme (b) Upwind moment scheme 


Fig. 8.3. Advection problem; 125 time steps with a = 0.8; distance travelled: 1 period; (a) solution by 
standard upwind scheme with Van Albada limiter; (b) solution by the moment scheme with the new 
limiter. 


The next test has both shocks and local extrema (Shu and Osher 1989). On [—1, 1], a moving Mach 3 
shock interacts with sine waves in density described by 

( (3.857,2.629,10.333), -1 < x < -0.8 

(p,u,p) = j . (7.9) 

( (1 + 0.2 sin(5n: x) , 0, 1), —0.8 < x < 1 

The final time is t = 0.36. Fig. 8.4 shows the results with 300 cells by standard upwind (with Van 
Albada limiter) and the upwind moment scheme. Here, for upwind moment, limiting is applied to the 
characteristic variables. The solid lines represent the solution with 3600 cells by the standard upwind 
scheme. Again, the upwind moment solution is highly accurate, and accuracy is preserved near extrema. 
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Fig. 8.4. Shu and Osher’s problem. The upwind moment solution is highly accurate: the shock is well 
captured and there is no loss of accuracy near extrema. 


The last test is the (2D) oblique shock problem (Yee, Warming, and Harten 1983). The domain is 
[0,4] X [0,1], The boundary conditions are: left, inflow with (p, u, v, p) = (l.,2.9, 0., 1/1.4); top, 
(1.7,2.6193, —.5063,1.5282); bottom, solid wall; right, outflow. The rectangular uniform mesh is 
80 X 20. Here, the limiter is applied in the x and y directions when needed (to the characteristic 
variables). Figure 8.5 shows (a) pressure distribution result by the upwind moment scheme and (b) 
pressures along y = 0.475 by the upwind moment and the standard upwind schemes. The solid line 
represents the exact solution. Again, the upwind moment solution is highly accurate. It appears that the 
slope update via moment calculation can ‘turn’ faster due to the compact stencil compared to the slope 
update via standard interpolation using neighboring data. For this case, compared with the standard 
(simplified) second-order upwind scheme with Van Albada’s limiter, computing cost increases by 
roughly a factor of five. 



Fig.8.5 (a) Oblique shock; pressure distribution by upwind moment scheme. 
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Fig.8.5 (b) Oblique shock solutions along y = 0.475 


4. Conclusions and Discussion 

In conclusion, two new high-order space-time schemes for conservation laws were introduced: the first 
is explicit and the second, implicit. The explicit method called the moment scheme achieves a CFL 
condition of 1 for the case of one-spatial dimension regardless of the degree of approximation. Fourier 
analyses for the one and two-dimensional cases were carried out. The moment scheme is accurate to order 
2p + 1 if the solution is approximated by piecewise polynomials of degree p. The implicit method is a 
simplified but optimal version of the discontinuous Galerkin scheme applied to time discretization. It 
reduces to a collocation implicit Runge-Kutta method using the right Radau points called Radau 11A. The 
explicit and implicit schemes are closely related since they employ the same intermediate time levels, and 
the former can serve as a key building block in an iterative procedure for the latter. A limiting technique 
for the piecewise linear scheme that can suppress oscillations near discontinuities while preserving 
accuracy near extrema was presented. Preliminary numerical results were shown. 

Whereas the current findings are encouraging, further research on both time-stepping and limiting is 
needed. 
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